
#include <stdlib.h>
#include <time.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <stdio.h>
#include <gl\glew.h>
#include <gl\freeglut.h>
#include <iostream>
#include <vector>
#include "kernel.cuh"




#define MAX_X      500
#define MAX_Y      500
#define SIZE	   MAX_X * MAX_Y
#define DIFF       0.0
#define DELTA_TIME 0.1

#define IX(i,j) ((i * MAX_Y) + (j))

using namespace std;

struct Vertex
{
    GLfloat x,y,red,green,blue;
};

Vertex *arr_h = 0;
Vertex *arr_d = 0;

float *density_d     = 0;
float *densityPrev_d = 0;
float *velocityX_d   = 0;
float *velocityY_d   = 0;
float *velocityX_h   = 0;
float *velocityY_h   = 0;

const size_t SIZE_BYTES = sizeof(Vertex) * SIZE;

dim3 blocksPerGrid(MAX_X, MAX_Y);
dim3 threadsPerBlock(1);
bool click;
int clicked_x;
int clicked_y;

unsigned int count = 0;

__global__ void kernel_part1(float *density_d, float *densityPrev_d, float delta_time, int size_x, float * velocityX_d, float * velocityY_d, Vertex * arr_d, int dir)
{
	
	int i = (blockIdx.x * blockDim.x) + threadIdx.x;
	int j = (blockIdx.y * blockDim.y) + threadIdx.y;
	int index = IX(i,j);

	density_d[index] -= delta_time * densityPrev_d[index];
	density_d[index] = density_d[index] * 0.69;

	//density_d[index] += delta_time * densityPrev_d[index]; //add fields

	float temp  = density_d[index]; //swap
	density_d[index]  = densityPrev_d[index];
	densityPrev_d[index] = temp;


	density_d[index] = densityPrev_d[index];//""linear solve"" 
	//advect
	

	//float x,y, s0, t0, s1, t1, dt0;

	//dt0 = delta_time * size_x;

	//x = i - dt0 * velocityX_d[index];
	//y = j - dt0 * velocityY_d[index];
	//if (x < 0.5f)        { x = 0.5f; }
	//if (x > size_x + 0.5f) { x = size_x + 0.5f; }
	//i0 = (int)x; 
	//i1 = i0 + 1;
	//if (y < 0.5f)        
	//{
	//	y = 0.5f; 
	//}
	//if (y > size_x + 0.5f) 
	//{
	//	y = size_x + 0.5f; 
	//}
	//j0 = (int)y;
	//j1 = j0 + 1;

	//s1 = 0;//x - i0;
	//s0 = 1 - s1;
	//t1 = 0;//y - j0;
	//t0 = 1 - t1;
	//		
	//density_d[index] = s0 * (t0 * densityPrev_d[IX(i0,j0)]);// + t1 * densityPrev_d[IX(i0,j1)]);// + s1 * (t0 * densityPrev_d[IX(i1,j0)] + t1 * densityPrev_d[IX(i1,j1)]);
	float dt;	
	int backjected_x, backjected_y;

	dt = delta_time * size_x;

	
	backjected_x = i - dt * 0.5 * velocityX_d[index];
	backjected_y = j - dt * 0.5 * velocityY_d[index];
			
			
	
	if (backjected_x < 0.5f)        
		backjected_x = 0.5f;
	if (backjected_x > size_x + 0.5f) 
		backjected_x = size_x + 0.5f;
		
	if (backjected_y < 0.5f)        
		backjected_y = 0.5f; 
	if (backjected_y > size_x + 0.5f) 
		backjected_y = size_x + 0.5f; 
	
	////s1 = new_pos_x - new_pos_x;
	//s1= 0;
	//s0 = 1 ;//- s1;
	////t1 = new_pos_y - new_pos_y;
	//t1 =  0;
	//t0 = 1;// - t1;	
	//density_d[index] = s0 * (t0 * densityPrev_d[IX(new_pos_x,new_pos_y)]);// + t1 * densityPrev_d[IX(new_pos_x , new_pos_y+1)]) + s1 * (t0 * densityPrev_d[IX(new_pos_x + 1 , new_pos_y)] + t1 * densityPrev_d[IX(new_pos_x + 1 , new_pos_y + 1)]);
	
	//density_d[index] = densityPrev_d[IX(backjected_x,backjected_y)] + densityPrev_d[IX(backjected_x , backjected_y+1)] +densityPrev_d[IX(backjected_x + 1 , backjected_y)] + densityPrev_d[IX(backjected_x + 1 , backjected_y + 1)];
	density_d[index] =  (densityPrev_d[IX(backjected_x,backjected_y)] + densityPrev_d[IX(backjected_x , backjected_y+1)] + densityPrev_d[IX(backjected_x + 1 , backjected_y)] + densityPrev_d[IX(backjected_x + 1 , backjected_y + 1)]) ;

	//coloring time
	float inten = density_d[index];

	//arr_d[i*MAX_Y + j].red   = inten/4;
	//arr_d[i*MAX_Y + j].blue  = inten/4*3;
	//arr_d[i*MAX_Y + j].green = 0;

	
	arr_d[i + j*MAX_Y].red   = inten/4;
	arr_d[i + j*MAX_Y].blue  = inten/4*3;
	arr_d[i + j*MAX_Y].green = 0;

}


__global__ void initArray(float *array, float value ) {
	int i = (blockIdx.x * blockDim.x) + threadIdx.x;
	int j = (blockIdx.y * blockDim.y) + threadIdx.y;

	//array[i*MAX_Y + j] = value;
	array[i + j*MAX_Y] = value;
}

__global__ void initArray2(float *array, float value, int point_x, int point_y ) {
	int i = (blockIdx.x * blockDim.x) + threadIdx.x;
	int j = (blockIdx.y * blockDim.y) + threadIdx.y;

	if (i > point_x && j > point_y && i < point_x+2 && j < point_y+2) {
		//array[i*MAX_Y + j] = value;
		array[i + j*MAX_Y] = value;
	}
}
__global__ void initArrayVelocities(float *array)
{
	int i = (blockIdx.x * blockDim.x) + threadIdx.x;
	int j = (blockIdx.y * blockDim.y) + threadIdx.y;
	array[i*MAX_Y + j] = 0.5;

}

void calculateDensity(int x_clicked, int y_clicked, int dir) {
	static unsigned int frame = 0 ;
	frame++;

	//if (!(frame%30))
	if (frame == 1)
	{
		initArray2<<<blocksPerGrid, threadsPerBlock>>>(densityPrev_d, 0.5,10,10);
	}
	else if(click == true)
	{
		printf("%d %d\n",x_clicked, y_clicked);
		click = false;
		initArray2<<<blocksPerGrid, threadsPerBlock>>>(densityPrev_d, 0.5,x_clicked,y_clicked);
		
	}
	
	//Sleep(800);
	//Sleep(750);
	kernel_part1<<<blocksPerGrid, threadsPerBlock>>>(density_d, densityPrev_d, DELTA_TIME, MAX_X, velocityX_d, velocityY_d, arr_d, dir);


	cudaError_t cudaStatus = cudaDeviceSynchronize();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaDeviceSynchronize returned error code %d after kernel!\n", cudaStatus);
    }

	cudaStatus = cudaMemcpy(arr_h, arr_d, SIZE_BYTES, cudaMemcpyDeviceToHost);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMemcpy failed!");
    }

}

int main_cuda(int argc, char **argv)
{
	arr_h = (Vertex *)malloc(SIZE_BYTES);
	for (int i = 0 ; i < MAX_X ; i++)
	{
			for(int j = 0; j<MAX_Y;j++)
			{
				arr_h[i*MAX_Y + j].x     = i;
				arr_h[i*MAX_Y + j].y     = j;
				arr_h[i*MAX_Y + j].blue  = 0;
				arr_h[i*MAX_Y + j].red   = 0;
				arr_h[i*MAX_Y + j].green = 0;

			}
	}

	// allocate array on device 
	cudaError_t cudaStatus = cudaMalloc((void **) &arr_d, SIZE_BYTES);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
    }

	// Density_d
	cudaStatus = cudaMalloc((void **) &density_d, sizeof(float) * SIZE);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
    }
	
	initArray<<<blocksPerGrid, threadsPerBlock>>>(density_d, 0.0);
	//initArray2<<<blocksPerGrid, threadsPerBlock>>>(density_d, 1);

	// DensityPrev_d
	cudaStatus = cudaMalloc((void **) &densityPrev_d, sizeof(float) * SIZE);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
    }

	initArray<<<blocksPerGrid, threadsPerBlock>>>(densityPrev_d, 0);
	//initArray2<<<blocksPerGrid, threadsPerBlock>>>(densityPrev_d, 0.5);

	// VelocityX_d
	cudaStatus = cudaMalloc((void **) &velocityX_d, sizeof(float) * SIZE);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
    }

	srand(time(NULL));


	 velocityX_h = (float *)malloc(SIZE * sizeof(float));
	 velocityY_h = (float *)malloc(SIZE * sizeof(float));
	 for (int i = 0; i < (SIZE);i++)
	 {
		 velocityX_h[i] = ((float)rand()/(float)RAND_MAX)/5;
		 velocityY_h[i] = ((float)rand()/(float)RAND_MAX)/10;
	 }


	
	initArrayVelocities<<<blocksPerGrid, threadsPerBlock>>>(velocityX_d);

	// VelocityY_d
	cudaStatus = cudaMalloc((void **) &velocityY_d, sizeof(float) * SIZE);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMalloc failed!");
    }
	
	initArrayVelocities<<<blocksPerGrid, threadsPerBlock>>>(velocityY_d);


	cudaStatus = cudaMemcpy(arr_d, arr_h, SIZE_BYTES, cudaMemcpyHostToDevice);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMemcpy failed!");
    }

		 cudaStatus = cudaMemcpy(velocityX_d, velocityX_h, SIZE * sizeof(float), cudaMemcpyHostToDevice);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMemcpy failed!");
    }

	cudaStatus = cudaMemcpy(velocityY_d, velocityY_h, SIZE * sizeof(float), cudaMemcpyHostToDevice);
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaMemcpy failed!");
    }
    // cudaDeviceReset must be called before exiting in order for profiling and
    // tracing tools such as Nsight and Visual Profiler to show complete traces.
    /*cudaStatus = cudaDeviceReset();
    if (cudaStatus != cudaSuccess) {
        fprintf(stderr, "cudaDeviceReset failed!");
        return 1;
    }*/
	
    return 0;
}